#### Report various yield statistics from data and model fits

## Setup environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- T
source("load.R", encoding="iso-8859-1")

library(dplyr)

## Report results from dataset

muni.arabica <- unique(df.variety$Munic�pio[df.variety$arabica.portion == 1])
muni.robusta <- unique(df.variety$Munic�pio[df.variety$arabica.portion == 0])

ii.arabica <- which(df$Munic�pio %in% muni.arabica)
ii.robusta <- which(df$Munic�pio %in% muni.robusta)

obsyield.arabica <- mean(df$total.produce[ii.arabica] / df$total.harvest[ii.arabica], na.rm=T)
obsyield.robusta <- mean(df$total.produce[ii.robusta] / df$total.harvest[ii.robusta], na.rm=T)

## Report results from model fits

sumstats <- read.csv("sumstats-eq24-orig.csv")
names(sumstats)[1] <- "fullparam"
sumstats$param <- gsub("\\[\\d+\\]", "", gsub("_coeff\\[(\\d+)\\]", "_coeff\\1", sumstats$fullparam))

truyield.arabica <- mean(sumstats$mean[sumstats$region %in% muni.arabica & sumstats$param == 'yield_mean'])
truyield.robusta <- mean(sumstats$mean[sumstats$region %in% muni.robusta & sumstats$param == 'yield_mean'])

truyield.arabica / obsyield.arabica
truyield.robusta / obsyield.robusta

death.arabica <- mean(sumstats$mean[sumstats$region %in% muni.arabica & sumstats$param == 'death_part'])
death.robusta <- mean(sumstats$mean[sumstats$region %in% muni.robusta & sumstats$param == 'death_part'])

range(df$year[!is.na(df$edd1000)])
avgprice <- colMeans(prices[prices$year >= 1980 & prices$year < 2018,])

truyield.arabica * avgprice[2]
truyield.robusta * avgprice[3]
